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Abstract 


A general strategy exists for constructing Energy stable Weighted Essentially TVon- 
Ctecillatory (ESWENO) finite difference schemes up to eighth-order on periodic do- 
mains. These ESWENO schemes satisfy an energy norm stability proof for both 
continuous and discontinuous solutions of systems of linear hyperbolic equations. 
Herein, boundary closures are developed for the fourth-order ESWENO scheme that 
maintain wherever possible the WENO stencil biasing properties, while satisfying 
the summation-by-parts (SBP) operator convention, thereby ensuring stability in 
an L 2 norm. Second-order, and third-order boundary closures are developed that 
achieve stability in diagonal and block norms, respectively. The global accuracy 
for the second-order closures is three, and for the third-order closures is four. A 
novel set of non-uniform flux interpolation points is necessary near the boundaries 
to simultaneously achieve 1) accuracy, 2) the SBP convention, and 3) WENO stencil 
biasing mechanics. 
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1 Introduction 

High order Weighted Essentially IVon- Oscillatory (WENO) methods are ideally 
suited for high fidelity simulations of complex physics containing discontinuities. 
An example of a typical application is a canonical simulation of sound generated by 
a shock-vortex interaction [1,2]. The high-order nature of WENO efficiently resolves 
the subtle details of the sound generation and propagation, while the stencil biasing 
mechanics ensures robust, high resolution properties in the vicinity of the shock. 

Although ideal for infinite or periodic domains, conventional WENO formula- 
tions encounter serious challenges at nodes that are near domain boundaries. To 
illustrate this, first realize that all high-order finite-difference (HOFD) formulations 
use “inward biased” stencils near boundaries that maintain the accuracy and sta- 
bility of the interior scheme. Next, note that WENO schemes invoke stencil biasing 
mechanics throughout the entire domain, including nodes adjoining the boundaries 
(e.g. fourth-order WENO schemes test three candidate stencils at each node). Thus, 
not only must the boundary closures used for WENO schemes be 1) inward-biased, 
2) stable and 3) accurate, they must also be stable for any possible combination 
of candidate stencils occurring anywhere in the domain. Simultaneously satisfying 
these constraints is a remote possibility if left strictly to chance. 

Many ad hoc boundary closures have been proposed for conventional WENO 
schemes, but most implementations still rely on low order boundary closures to 
achieve stability and robustness. For flow near a wall, Shen et al. [3] constructed the 
domain such that the no slip condition is enforced at the WENO flux point. Pressure 
was extrapolated to the wall using a third order scheme, and the next WENO flux 
cell was treated with a third order reconstruction. This procedure reduces the global 
order of accuracy of their 5th order WENO scheme. Alternatively, ghost points can 
be used in combination with physical boundary conditions to specify data in stencils 
that extend outside of the domain [4] . This reduces the generality of the boundary 
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treatment. Another strategy is to simply discard stencils that extend outside of the 
domain [4]. 

The principle objective of this work is to develop boundary closures to com- 
plement the periodic domain Energy .Stable WENO finite difference methodology 
developed and reported in Refs. [5] and [6]. The new interior/boundary ESWENO 
schemes (henceforth referred to as “finite- domain ESWENO”) retain all the salient 
features of the original periodic schemes, including: 1) conservation and L 2 -energy 
stability for constant coefficient (linear) hyperbolic systems, including those with 
discontinuous initial or boundary data, 2) design order accuracy throughout the 
entire domain, especially regions near the boundaries or near smooth extrema, and 
3) full stencil biasing mechanics (Left, Central, Right) at all possible points. 

The finite-domain ESWENO schemes are constructed by adding a “special” non- 
linear artificial dissipation term to a conventional WENO scheme (e.g. those found 
in Refs. [7,8]). The additional term is design-order accurate for smooth solutions, 
including smooth extrema, and is constructed such that the resulting ESWENO 
scheme is stable in the L 2 -energy norm for both continuous and discontinuous so- 
lutions. 1 In this initial work, we focus exclusively on the central fourth-order class 
of ESWENO schemes, leaving the development of 6th- and 8th-order schemes to a 
later date. 

The strategy for developing finite-domain ESWENO schemes is for the most 
part equivalent to that used to build the periodic domain ESWENO schemes. The 
essential elements are summarized as follows: 

• Develop a 4th-order finite-domain target scheme that is stable, conservative 
and accurate for smooth flows. This task is accomplished using the summation- 
by-parts (SBP) [6,9] framework (as was done in the periodic case). The SBP 
framework requires the construction of accurate discrete schemes that satisfy 
basic matrix conditions, that if met guarantee discrete stability and conserva- 
tion. 

• The target scheme is recast into the “dual-grid” framework of the conventional 
WENO approach, whereby the solution is stored/advanced at the gridpoints, 
while interface fluxes are constructed at the “half-points” . 

• The stability of the resulting stencil is tested, and if unstable, is modified with 
a special artificial dissipation term. This ensures that the discrete eigenvalues 
of the spatial scheme are bounded within the left-half of the complex plane. 

The organization of the paper is as follows. Section 2 introduces a dual-grid 
nomenclature and summarizes the necessary SBP formulas including a discrete sta- 
bility proof. Section 3 summarizes the implementation of the ESWENO scheme, 
and then expresses the baseline ESWENO operator in SBP form. Section 4 formu- 
lates the artificial dissipation term necessary for the ESWENO approach. Section 

1 The conventional WENO scheme, does not have such an energy estimate. Thus, while the 
spectrum of the symmetric part of the ESWENO scheme is always located in the left half of the 
complex plane, the symmetric part of the WENO scheme could have positive eigenvalues. 
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5 presents a summary of numerical test cases. Conclusions and directions are pre- 
sented in Section 6. An appendix is included with all the implementation details. 


2 Definitions 

2.1 Complementary Grids 

The implementation of the ESWENO (WENO) stencil biasing mechanics, as a mat- 
ter of convenience, uses two complementary sets of gridpoints, differing in dimension 
by one. As such, two discrete sets of points 

X7V = [xi,X2,’ mm ,XN-1 ,Xn] T ; X M = [x 0 , X \ , • • • , XjV_i , X N ] T 

are defined on the interval 0 < x < 1. Furthermore, a projection is defined of any 
function (j) onto Xjv and 5cm respectively as 

4> = [<j)(xi ), 0(z 2 ), ■ • • , <j)(x N -l), 0 {xn)} T ; <i> = [<j>(x 0 ),<f)(x l), • ■ ■ , (/>(xjv-l), ^>(xat)] T 


The xtv is a uniform set of “solution points” and stores the discrete solution u. 
Surrounding each solution point, xjvj, a control volume is constructed, whereby the 
governing conservation laws are enforced. A second set of “flux points”, xjf. is sit- 
uated at the interfaces between adjoining control volumes, where data interpolated 
from u is used to construct interface fluxes / and establish a simple flux balance. 

Figure (1) illustrates the relationship between the solution-nresh x^r and the 
flux-mesh 5cm- 

r r 


fi—l fi fi + 1 fi+2 



Sl S c S r 


Figure 1. The computational domain is depicted. The ESWENO solution and flux 
points, and the smoothness indicator stencils S L , S c and S R are shown. 

The biasing mechanics move data via interpolation or differentiation between 
the two sets of points using non-square matrix operations. For example, a flux 
constructed on the high- dimensional set from data on the low-dimensional set, is 
represented as 

f = m%n f (1) 

where M = N + 1 and m^n is an [(A" + 1) x N] matrix. 
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2.2 Summation-By-Parts Operators 

Consider a linear, scalar wave equation 

If + m = °> f = au , t > 0, 0 < x < 1, ^ 

rt(x,0) = rt 0 (x) ; u(0,f) = g(i) 

where u is the continuous solution, a is a constant and uo(x) is a bounded piecewise 
continuous function. Without loss of generality, assume that a > 0 on the interval 
0 < x < 1. Applying the energy method to equation (2) leads to 

JdMlL + a[u 2 {l,t) ~g 2 (t)] =0 (3) 

where || ■ \\l 2 is the continuous L 2 norm. Thus, equation (2) admits an energy solution 
that is strictly dissipative (if g(t) = 0, a > 0) for all time. We now develop using 
mimetic techniques (see Ref. [6] or Ref. [10]) a class of discrete spatial operators for 
which the discrete energy is bounded from above by the continuous target estimate 
provided by equations (2) and (3). 

Define the low dimensional space x^v to be a uniform grid Xj+i = jAx, j = 
0, • • ■ , N — 1, with Ax = 1 /{N — 1). On this grid, we define a flux f = au and 
its derivative f x = au x , where u = [u(xi, t), . . . , u(x at , t)] T and u x = [u x (x\, t), . . . , 
Ux(xn, t)] T are projections of the continuous solution u and its derivative u x onto the 
computational grid x^r- Next, define a fourth-order approximation for the first-order 
derivative term in equation (2) as 

— = Df + 0(Ax 4 ) (4) 

Placing a mild restriction on the generality of the derivative operator (see Ref. [10] 
or Ref. [6]), the matrix D is expressed in the following form: 


D = V~ 1 [Q + K) ; Q + Q T = Diag[- 1,0, ■■■,0,1] 

n = lZ T ; v t Kv> 0 (5) 

V = V T ; v T Pv > 0 

for any real vector v/0. The matrix TZ is defined by the expression 

n = A 0 + D\Ai[Di] t + D\ [Di] t A 2 D 1 [Di\ t + D 1 [D 1 ] T D 1 A 3 [D 1 ] T D 1 [D 1 ] r (6) 

which allows a constructive means of forming a symmetric semi-definite discrete 
operator. Herein we define the matrices A* to be diagonal positive semidefinite 
matrices of the appropriate size: i.e. 

A e = Diag[Xi, ■ ■ ■ , X N ] ; Xj > 0 , j = 1, N ; e = 2p,p = 0, 1 

A 0 = Diag[\ 0 , ■ ■ ■ , Xn] ; Xj>0,j = 0,N ; o = 2p+l,p = 0 , 1 
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Using the definitions of the SBP operators [eqs. (4) — (6)] , and defining on x^r 
the discrete approximate solution u and flux f = au, the semi-discrete counterpart 
of equation (2) becomes 


-V 


-l 


9u 

dt 


p- 


Q f = 

T 


A 0 + D 1 A 1 [D 1 ] t + D^Drf A 2 D 1 [D 1 ] 


U 1 [U 1 ] t ^ 1 A 3 [U 1 ] t ^ 1 [D 1 ] t 


(7) 


To show that the above finite-difference scheme is stable, the energy method is 
used. Multiplying equation (7) with u T P, adding the result to its transpose, and 
rearranging the expression using equation (5), yields the semi-discrete energy equa- 
tion 2 


4ll u llp + a[u 2 (N,t) -g 2 (t)\ = -2a^2([D 

i = 0 


iiT 


u) T A i[D 


ii T u < 0 


( 8 ) 


The right-hand side of equation (8) is non-positive because the diagonal matrix 
A n is positive semidefinite [v T A n v > 0 for all real v of length N] and a > 0; thus 
the stability of the finite-difference scheme given by equation (7) is assured. This 
result can be summarized in the following theorem: 


Theorem 1 The approximation [eq. (7)/ of the problem [eq. {2)] is stable if equa- 
tions / (4) —(6) 7 hold. 


Remark. The finite-difference scheme [eq. (7)] constructed as an approximation 
of equation (2) is inherently nonlinear, despite the linearity of (2). Indeed, the 
matrices Q and A n are nonlinear matrices [i.e. , Q = Q(u) and A n = A n (u)]. 

Remark. Although the matrices Q and A n are constrained in form by Q + Q T = 
Diag[— 1,0, ■ ■ • , 0, 1] and v T A n v > 0, no specific assumptions are made about the 
values of the coefficients. 

Remark. All matrices present in equation (7) need not be square. Indeed, it is 
advantageous to define the [N x (IV + 1)] dimensional matrix D \ to be 


D i = 


/ "I 
0 

0 
0 
0 


1 

-1 

0 

0 

0 


0 

1 

0 

0 


0 0 0 \ 

0 0 0 

0 0 

-1 1 0 

0 — 11 / 


(9) 


to properly account for the difference between the dimension of the x^r and xjvf 
grids. 


2.2.1 Central Fourth-order SBP Operators 

The target operator used herein for the ESWENO scheme is a fourth-order non- 
dissipative (i.e. 1Z = 0) central SBP operator, for which equation (5) reduces 

2 The nomenclature Dd, i = 1, q used in equation (8) is not valid for nonsquare matrices. Herein, 

define D\ 2p to mean {Di[Di] t } 7 ; p = 1, q. 
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to D = V X Q subject to the constraints Q + Q 1 = Diag[— 1,0, ••• ,0,1] and 
V = V T ; v T Pv > 0. 

The matrices V and Q are composed of the conventional banded operators in 
the interior of the domain: 


Pj = Ax I 


Qi 


Pentadiagonal 


1-8 8 - 1 ' 

12 ’ 12 ’ 12 ’ 12 


Boundary closures satisfying the SBP convention require auxiliary stencils at Nb 
points near each boundary. Elsewhere, it is shown [11] that iVj, > 4 is required to 
satisfy the necessary accuracy, symmetry and skew-symmetry constraints. Herein, 
the Nb = 4 is assumed, whereby the matrices V and Q take the form 


( p o 

0 

0 \ 

1 

( Qo 

Qd 

0 \ 

O 

B 

< 

II 

I 

0 

; Q = 

-Q T d 

Qi 

Qd 

V ° 

0 

(Po) PT ) 

\ 

\ o 

~Qd 

- m PT ) 


(10) 


with 


Po = 


pn 

P 1 2 

P13 

PlA 

\ 


( 

-1 

2 

712 

713 

714 

\ 

P12 

P22 

P23 

P24 


; Qo = 


Nl2 

0 

723 

724 


P13 

P23 

P33 

P34 




-<713 

“723 

0 

734 


Pl4 

P24 

P34 

PU 



v - 

-714 

— 724 

“734 

0 

/ 








f 0 












o 










Qd 


-1 

n 











12 

u 











\ — 
\ 12 

-l 

12 

0 • 

•• ) 



( 11 ) 


where PT denotes the per-symmetric transpose. 

Two types of P-norms are typically used to prove stability with SBP discretiza- 
tions: diagonal, and block. The diagonal-norm is advantageous when considering 
equations of more complexity than the constant coefficient, linear equation con- 
sidered herein. Indeed, SBP proofs using the diagonal-norm extend to variable- 
coefficient problems [12], as well as certain nonlinear formulations [13,14], a property 
not shared by the block-norm boundary closures. Diagonal-norm boundary closures 
can not preserve the full design-order potential of the interior operator. Closures 
satisfying the matrix constraints in equation (8) with a fourth-order interior scheme, 
can have accuracy no higher than second-order. Thus, the global accuracy of this 
combination is only third-order [15, 16], or in general p + 1-order despite 2p-order 
accuracy in the interior. (The nomenclature p — 2p — p\ p = 1, • • • is sometimes used 
to describe the diagonal-norm formulations. See Ref. [9] and Ref. [11] for further 
discussion.) 

Although block-norm closures are less robust than their diagonal-norm coun- 
terparts, for variable-coefficient linear or nonlinear equations, they enable the use 
of 2 p — 1-order stencils, and preserve the formal accuracy of the 2p-order interior 
operator. For example, third-order boundary closures preserve the formal accuracy 
of the fourth-order interior scheme. Furthermore, the additional parameters in the 
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■P-norm (six, in the (4 x 4) case) can be used to optimize other properties of the 
third-order closure. Indeed, two free parameters (herein denoted i9/, l = 1,2) remain, 
after satisfying all accuracy, and SBP matrix constraints, when using a (4 x 4) block 
norm to close a fourth-order interior. Both parameters, as will be shown in the 
next section, are needed to develop stable, and accurate fourth-order ESWENO 
boundary closures. 


3 Expressing Summation-By-Parts Operators in Flux 
Form 

The next objective is to recast the target fourth-order central SBP operator ( D = 
V~ 1 Q defined on the x^r grid) into the “dual-grid” framework of the conventional 
ESWENO/WENO approach. Of particular interest is extending the ESWENO’s 
“interior” stencil biasing mechanics, to those stencils at or near domain boundaries. 
We begin with a cursory description of the fourth-order ESWENO scheme, followed 
by its representation in matrix form. 


3.1 General Structure of ESWENO Schemes 


The conventional high-order WENO/ESWENO finite difference scheme for the scalar 
1-D wave equation (Eq. (2)) can be written component-wise in the following semi- 
discrete form: 


duj , 4+§ 4'-§ n 

—r~ H wz = 0 

at AxMj 


(12) 


The numerical flux /• , i for the central fourth-order WENO scheme is given by 
the expression 


4+5 W j+l/2fj+l/2 + W f+l/2fj+l/2 + W f+l/2fj+l/2’ ( 13 ) 


where f( r ) = m4v h r = {h, C, R} are second order fluxes constructed from data 
interpolated using three distinct stencils S L , S c , and S R . 3 

The nonlinear weight functions j , 2 • r = {L,C,R} used in equation (13), 
are constructed such that the interface flux /■ , i is a convex combination of the 
three candidate second-order fluxes. The weight functions are constructed using the 
expressions 


w 


(r) 

3 + 1/2 


_(r) 

a j+ 1/2 

E -(r) 
a 3 + 1/2 


_ j{r) 

5+1/2 U j+ 1/2 



T J+l/2 \ 

6 + Pj r +l /2j 


r = L,C,R (14) 


^Figure ( 1 ) shows a schematic of the interpolation process. Each second-order flux is con- 
structed from data obtained using only two points: S L = {xj-i,Xj}, S c = {xj,Xj+ 1}, and 
S R = {xj+i, Xj+2} ■ The r = {L,C,R} nomenclature identifies the origin of the data relative to 
the interface position. The interpolates need not be limited to three; indeed S LL = 
and S RR = {xj+2, £5+3} can be included near boundaries to increase the fidelity of the boundary 
closures. 


~(r) 

where d\^ , 2 is the target weight of the candidate stencil. If no stencil biasing occurs 
(a well resolved solution), the target central difference operators are recovered, and 
the functions w ^ used as weights in equation (13) reduce to 


-M 
W j+ 1/2 


cf r) 

j+1/2 
J + 1/2 


d ir) 


j+1/2' 


- ~(r) 

Definitions of the parameters Tj +1 / 2 , e, /T + j , 2 , r = L, C, R, as well as implementation 
details are included in the appendix. 

The flux derivative expressed in equation (12) in component form can be written 
as 

' J 1 0 i - (pm-af), as) 

where [<5 xm] is the diagonal matrix [<5 xm] = Diag[DjXM] that accounts for variable 
cell sizes. Further expanding the flux function f in equation (15) in terms of matrix 
interpolation operators, and equating the result with the target SBP form, yields 
the expression 

= [65c m ]~ 1 D 1 [i Cu L (d ] ) m Tn + “ C (d) M ^ + uo R {d) m In\ f = P _1 Qf 

( 16 ) 

The solution to equation (16) is the central objective of this section; specifically, to 
express a target fourth-order SBP discretization V~ 1 Q, in terms of the conventional 
1) interpolation operators m^n, 2) weight functions Co, and 3) difference operator 
D\ needed to implement the ESWENO formulation. 

Remark. There is no guarantee that a solution exists to Equation (16). The 
conventional periodic domain operators m^ t n , Co r ; r = L,C, R on uniformly spaced 
grids xtv , xm, match the fourth-order centered 'P~ 1 Q identically in the interior of 
the domain. The task at hand is to determine boundary operators for n , Co r ;r = 
L,C,R and the grids x^r , xjf that match T~ l Q at boundaries. 

Remark. Equation (16) is a nonlinear function of the flux grid Indeed, 

assuming a uniform solution grid x/y, the functional dependency of the matrix op- 
erators are 



dW(xAf); M lP = a 4 r) (* M ) 




r = L,C,R 


T = V(0i), Q = Q($/) , 1 = 1,2 


3.2 Nonuniform Flux-Points 

An analytic solution to equation (16) is not immediately forthcoming using Mathe- 
rnatica 7.0 [17], if the solution xjy and flux xm points are uniformly spaced on the 
domain 0 < x < 1: 


x j+ 1 
Xk 


jAx , j = 0, N — 1 , Ax = 1/(N — 1) ; 

o x k +x k+1 ) t k = l,N- 1, x 0 = 0,x N = l, 
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for either the diagonal-norm or the block-norm SBP operators V~ 1 Q. Increasing the 
boundary block sizes to (6x6), while providing numerous additional free parameters 
di, does not remove an apparent inconsistency in the constraint equations. 

Multiple solutions exist for equation (16), for either the diagonal- or block-norms, 
if the second, third, and fourth flux points xm are shifted to a nonuniform set of 
points near the boundaries. (Figure (1) illustrates the relationship between the 
solution-mesh x^r and flux-mesh xjf, including the shifted flux points near each 
boundary.) The location of the flux points is precisely related to the P-norm used to 
prove stability for the SBP operator. The following lemma is a central contribution 
of this work: 

Lemma 1 Assume N uniformly distributed solution points x^r . There exists an N + 
1 dimensional set of flux points xm for which the central , fourth-order, ESWENO 
interpolation operators, can be expressed in SBP form. The discrete set of points 
xm satisfies the following expression: 

V~ x D x *m = 1 (17) 


with 1 the unit vector. 


Proof: The proof is constructive. Given xm ; the general solutions are determined 
using Mathematica. Owing to the complexity of the general expressions, only two 
optimized solutions are included herein: one derived using a diagonal-norm, the 
other a block-norm. Note that D\5 cm = VI is an equivalent statement because V 
is nonsingular. □ 

Although a mathematical motivation for Lemma (1) has not been found, a di- 
mensional argument can be used to explain the relationship between V and xm- 
First, note that the D\ matrix is an undivided difference, and does not contain 
any reference to the precise location of the flux-points f, nor to their separation; 
D\ is dimensionless and simply forms the difference between the left and right cell 
interface states. The dimensionality (or inverse weight) in the derivative opera- 
tion is provided by the term [5xm] _1 - Left multiplying the flux difference D\f 
by 1 t [<1xm] (the matrix equivalent of integration), yields the domain flux balance 
l T [55tM\Dif = /( 1) — /( 0). Clearly, the individual terms in the matrix [(5 xm] are 
the quadrature weights needed for discrete integration. Thus, the matrix [<5xm] is 
filling precisely the same role in equation (16), as the V matrix in the SBP context 
D = V^ 1 Q. It is not surprising that V and [D\Xm\ are closely related. 

The next two examples demonstrate the validity of Lemma (1) as well as pro- 
vide the flux points xjf needed for the diagonal- and block-norm fourth-order SBP 
operators. First, consider the [N x N] diagonal P-norrn Vd defined by 

V, - A xDiaa\VL i? 43 « j ... j 1? « 5? lri 

/ d — i Jiuy 48 48 48 4 8 x x 48 48 48 48 J 


Distribute the flux points x^ as 


xm 


7Ax 

0,X±,X2,X 3 , — - 



1 T 


(1 - x 3 ) ,(1 -x 2 ) ,(1 - Xi) ,1 


(18) 
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where 


17Ax _ _ 59Ax _ _ 43Ax _ _ 49Ax 7Ax 

* 1 = * 0+ ^ 5 * 2 = * 1 + ^ 5 * 3 = * 2 + ^ 5 * 4 = * 3+ ^ = ^' 

Then, the vector D\5Lm is precisely the reciprocal of the diagonal element of Vd- 
The product V^ 1 D\5lm is trivially the identity vector 1, [the condition of Lemma 
(1)] since Vd is a diagonal matrix. 

Lemma (1) generalizes immediately to block-norm V cases. Recall that V is not 
uniquely defined in the fourth-order, centered, block-norm case. The matrix V is a 
function of two independent parameters (i.e. V = V(i ?;), Z = 1,2), after all accuracy, 
and matrix constraints are satisfied. In this specific case (but not in general) the 
matrix- vector product VI, used in Lemma (1) to determine xm, is independent of 
the parameters Thus, the following set of points xm 


(1 


XM = [0, 


43Ax 61Ax 349Ax 7Ax 


144 ’ 36 ’ 




144 > 2 

'#) , (1 


43Az\ I T 
144 / ’ 1 


(19) 


satisfied the conditions of Lemma (1), independent of the values of the available 
parameters $i, l = 1 , 2 . 

Remark. Existence of a set of points xm does not imply that xm is unique, 
although no other set has been identified. 

Remark. The solution- (xjv) and flux- (xm) grids are interdigitated, the excep- 
tions being the collocated endpoints. Given this set of flux points xm, the interpola- 
tion stencils S r , r = {L, C, R} are consistently defined for all points xmj ■ j = 2, N— 2. 
The flux-points next to each boundary xm ;] , j = 1 and j = N — 1 only have a centered 
and an inward-biased stencil. The first and last flux-points, in addition to inward 
biased interpolation stencils, have an “exact” interpolant derived from the collo- 
cated solution-point. This property is shared with conventional WENO boundary 
operators. 


3.3 Interpolation Operators and Preferred Weights 

An analytic solution to equation (16) exists for either the diagonal-norm or the 
block-norm case, when xm is given by equation (18) or (19), respectively. The 
diagonal-norm case required additional interpolants at the boundaries, to become 
realizable. Specifically, the interpolation operators I LL , I RR , and weights w LL , w RR , 
were added in equations (13) and (16). The additional terms were only required at 
the boundaries. 

Interpolation operators for the WENO/ESWENO 3 . 4.3 and WEN0/ESWEN0 2 _4 _ 2 
operators are given in the appendix. 


4 Energy Stabilization Terms 

When the nonlinear weights used in equation (13) satisfy the target conditions, 
= &( r \r = {L,C,R}, then all stencils are stable without the dissipation term 
V~ l 1Z. With nonlinear weights, however, there is no guarantee that the baseline 


11 


WENO operator D is stable. To ensure stability, an artificial dissipation term D a d is 
added to the WENO operator D, such that the symmetric portion of the resulting 
(ESWENO) operator D is positive-definite. Specifically, D satisfies the following 
matrix conditions: D = D + D a d subject to D sym + D a d > 0. 

The matrix D a d is constructed by expanding the symmetric matrix D sym into 
its elemental components using the decomposition 4 


D 


sym 


= v 


-i 


A 0 + D 1 A 1 [D 1 ] t + A 2 D l [D l ] T + . 


The existence of this decomposition is established in Ref. [6] . 

Next, the diagonal terms of Aj are modified to be smoothly positive using the 
expression: 


A?> = 1 


a?> r+i ?— a 


(0 


where 5 is dependent on the grid spacing. 

<5i = (A£) 3 6 2 = (AO 2 S 3 = (AO 2 

The artificial dissipation operator takes the final form 


( 20 ) 


(21) 


V 


-1 


Dadi = V-'Kf = 

Ao + D\Ai[Di] t + D l [D l ] T A 2 D 1 [D 1 ] T + 


f 

(22) 

The dissipation matrix Ao (at least herein) is identically zero for both the diag- 
onal and block norm schemes; thus it is convenient to define the vector 


4> = A^Thf + [E ) 1 ] T A 2 D 1 [I? 1 ] t + [Di] t D 1 A 3 [D\] T Di[Di] J 


(23) 


This definition of ^ allows the artificial dissipation term to be combined with the 
original flux f given by equation (13) to yield the following: 

^ = v- 1 D 1 (f + $) + 0(Ax 4 ) (24) 


5 Numerical Tests 


The accuracy and stability characteristics of the new ESWENO framework are tested 
using the one dimensional advection equation and the quasi-one dimensional and 
two dimensional Euler equations. As a more practical test, the framework was 
implemented in a two dimensional reacting Navier-Stokes solver and was used to 
simulate supersonic non-reacting and reacting shear layers. All simulations are 
integrated in time or pseudo-time using the same five-step, fourth-order Runge Kutta 
solver [18]. 

4 Note that the dimensions of the even indexed A j are (N x N), while the odd indexed A j have 
dimensions (N + 1 x IV + 1). This is the result of the difference matrix D 1 having dimensions 
[JVx (N + 1)]. 
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5.1 Eigenvalues of the ESWENO operator 


To demonstrate that the finite-domain ESWENO operator is positive semidefinite, 
an eigenvalue analysis is conducted using the one dimensional scalar advection equa- 
tion with discontinuous initial data. We begin by rearranging equation (2) into the 
form 

V d ^ = -a(Q + K) u (25) 

at 

Setting u = ue kt (see Ref. [19]) in equation (25) yields 

- Wu = (Q + K)u, A = - (26) 

a 

Stability is guaranteed if all real eigenvalues satisfy the condition A r < 0 since a > 0. 
ESWENO/WENO operators depend on the calculated solution, so eigenvalues were 
calculated over ten characteristic times using N = 100 grid cells to ensure consis- 
tency of the analysis. The boundary condition used in the simulation was g(t) = 0. 
Eigenvalues were computed for the ESWENO/WENOg,^ scheme applied to a 
convected square wave. The resulting maximum real eigenvalues are plotted against 
time in Figure (2). 



Figure 2. Maximum real eigenvalues for a convected square wave are plotted as a 
function of time for the WENOg^g operator with and without the energy stable 
addition. 

It is clear that the bounded WENC> 3 _ 4_3 operator has eigenvalues in the right 
half of the complex plane. This is fixed by the addition of the energy stable term in 
the ESWENC> 3 _ 4_3 scheme. This result is similar to that observed in the periodic 
case [5, 6]. 


13 


5.2 Linear Wave Equation 

Numerical solutions of the linear advection equation were examined to test the ac- 
curacy and dissipation characteristics of the finite-domain ESWENO schemes. The 
convergence rate of the ESWENC>2_4_2 and ESWENOg,^ operators are calculated 
for a smooth solution. To evaluate the robustness of the schemes, a wave solution 
with a time-dependent, discontinuous inlet condition is examined. It is important 
to note that since the ESWENOg,^ operator is constructed using a block P-norrn, 
simply assigning the boundary solution to be the boundary data (injection bound- 
ary conditions) does not maintain energy stability. Instead, both ESWENO schemes 
utilize the Simultaneous Approximation Term (SAT) penalty method to enforce the 
boundary conditions [20]. 

5.2.1 Sine Wave 

To test the global order of convergence, the ESWENO scheme is used to approximate 
the solution to the linear advection equation with the exact solution 

u(x,t ) = sin(x — at), a = 1, i£(- 7r,7r) (27) 

Ten characteristic times were simulated. The L 2 and L ^ errors are tabulated along 
with the convergence rate in Tables (1) and (2), respectively, for the full block-norm 
ESWENO operator and are compared to the errors from the block- norm linear cen- 
tered difference operator. For the diagonal norm operators, the errors are tabulated 
in Tables (3) and (4). It is clear that the ESWENO operator yields the design or- 
der of convergence for the linear advection equation, and provides almost the same 
accuracy as the linear counterpart. 

Table 1. L 2 error and convergence rate of a smooth solution to the ID advection 
equation for the linear and ESWENO operators with the block norm boundary 
closure , 



Linear Block Norm 

ESWEN0 3 _4_ 3 

Number of Cells 

L 2 Error 

L 2 Rate 

L 2 Error 

L 2 Rate 

100 

1.60e-06 

- 

1.65e-06 

- 

200 

1.01e-07 

3.99 

1.04e-07 

3.99 

400 

6.31e-09 

4.00 

6.55e-09 

3.99 

800 

3.95e-10 

4.00 

4.10e-10 

4.00 

1600 

2.52e-ll 

3.97 

2.61e-ll 

3.97 


5.2.2 Square Wave 

A square wave advection problem with the exact solution 

tanh (b (x — at + 2-7r)) — tanh (b (x — at + ^)) 
u{x, t) = — — ^ o=l, x £ ( — 7r, 7 r) 

(28) 
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Table 2. Loo error and convergence rate of a smooth solution to the ID advection 
equation for the linear and ESWENO operators with the block norm boundary 
closure, 



Linear Block Norm 

ESWENO 3 . 4.3 

Number of Cells 

Loo Error 

Loo ■R'Clt'G 

Loo Error 

Loo ILi-'tG 

100 

3.24e-06 

- 

3.43e-06 

- 

200 

2.07e-07 

3.97 

2.18e-07 

3.97 

400 

1.30e-08 

3.99 

1.38e-08 

3.99 

800 

8.17e-10 

3.99 

8.65e-10 

3.99 

1600 

5.20e-ll 

3.97 

5.51e-ll 

3.97 

Table 3. L 2 error and convergence rate of a smooth solution to the ID advection 
equation for the linear and ESWENO operators with the block norm boundary 
closure. 


Linear Diagonal Norm 

ESWEN0 2 _4_2 

Number of Cells 

L2 Error 

L2 Rate 

L2 Error 

L2 Rate 

100 

4.20e-05 

- 

2.05e-05 

- 

200 

5.21e-06 

3.01 

2.30e-06 

3.16 

400 

6.48e-07 

3.01 

2.72e-07 

3.08 

800 

8.09e-08 

3.00 

3.30e-08 

3.04 

1600 

1.01e-08 

3.00 

4.07e-09 

3.02 


was tested, where b = 200 caused a very sharp gradient that could not be resolved. 
The numerical solution was calculated for up to t = 10.0, and N = 200 grid cells were 
used. The initial condition was uq(x,0) = 0 in the domain. The wave was passed 
into the domain, advected through the interior, and then passed out of the domain. 
This problem tests the dissipation and oscillatory characteristics of the boundary 
treatment as well as the interior treatment. A comparison of the numerical solution 
to the exact solution is shown for the two different boundary norms in Figure (3) 
at two different times. The solutions shown in Figure (3) prevent oscillations, and 
smearing of the discontinuity caused by the boundary treatment, particularly at 
the inlet, is reasonable. The solution using the ESWENC> 2 _ 4_2 operator in the left 
figure (t = 2 . 0 ) exhibits slightly more dissipative error than the solution using the 
ESWENC> 3 _ 4_3 operator. 


5.3 Euler Equations 

The quasi-one dimensional Euler equations given by 



f p } 

, 1 

H 

pu 

K p e ) 

1 f =( 



G = 


pu 

puu 

{pE + p)u 


(29) 
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Table 4. Loo error and convergence rate of a smooth solution to the ID advection 
equation for the linear and ESWENO operators with the block norm boundary 
closure, 



Linear Diagonal Norm 

ESWEN0 2 _4_2 

Number of Cells 

Loo Error 

Loo ■R'Clt'G 

Loo Error 

Loo 

100 

7.79e-05 

- 

6.60e-05 

- 

200 

9.96e-06 

2.97 

8.16e-06 

3.02 

400 

1.25e-06 

3.00 

1.07e-06 

2.93 

800 

1.57e-07 

3.00 

1.37e-07 

2.97 

1600 

1.96e-08 

3.00 

1.73e-08 

2.99 



Figure 3. Numerical solutions using the ESWENC> 2 _ 4_2 and ESWENOg.^ oper- 
ator compared to the exact solution of an advected square wave at two different 
times. 


were used to test the extension of the ESWENO operators to non-linear hyperbolic 
systems of equations. In order to properly calculate stencil biasing parameters, 
the conservative fluxes in equation (29) are transformed into characteristic form 
and split into forward and backward propagating waves using global Lax-Friedrichs 
splitting. 

= - (F c ± A maa ,U c ) (30) 

where F c and U c represent the characteristic fluxes and variables, respectively. The 
transformation from conservative to characteristic variables is performed at each flux 
point using transformation matrices constructed from the Roe averaged variables at 
the flux point. The ESWENO reconstruction is used to calculate interpolated fluxes, 
and then the interpolated fluxes are transformed back to physical space, where the 
gradient is calculated. 

The solutions to the quasi- ID Euler equations utilize the SAT penalty boundary 
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conditions in order to ensure stability of the solution and consistency of the boundary 
closure. The details of the SAT penalty method for the Euler and Navier-Stokes 
equations can be found in Svard, Carpenter, and Nordstrom [21]. 

5.3.1 Subsonic One-Dimensional Nozzle 

The first nonlinear test problem was a steady state, subsonic flow in a one-dimensional 
nozzle. The area was defined as 

A{x) = 1 — 0.8x(l — x), xE(0, 1) (31) 

The inlet Mach number was set to M = 0.5 and the outlet pressure was set to 
be equivalent to the inlet. The exact solution is plotted in Figure (4). The L 2 
error norm of the density and the convergence rate for both ESWENO schemes is 
tabulated in Table (5). The calculated convergence rates match the design order. 

Remark. The design order of convergence is achieved despite the fact that the 
transformation matrices are calculated from second order accurate Roe averages. 
This is because the transformation is only used to transform to wave space and 
back to physical space in order to facilitate flux splitting. The accuracy of this 
transformation does not affect the order of accuracy of the interpolation. Instead, 
the accuracy of the transformation relates to how accurately the physical fluxes can 
be decomposed into waves. 

Table 5. L 2 error and convergence rate of a smooth solution to the quasi-one di- 
mensional Euler equations for ESWENC>2_4_2 and ESWENOg,^. 


Number of Cells 

ES WEN Oq_4_2 

ESWEN0 3 _4_ 3 

L 2 Error 

L 2 Rate 

L 2 Error 

L 2 Rate 

25 

6.55E-05 

- 

2.13E-06 

- 

50 

6.14E-06 

3.42 

1.14E-07 

4.22 

100 

8.05E-07 

2.93 

6.43E-09 

4.15 

200 

1.03E-07 

2.97 

3.82E-10 

4.07 


5.3.2 Transonic One-Dimensional Nozzle 

The next test of the boundary-closed ESWENO scheme was a steady state transonic 
one-dimensional nozzle. The cross sectional area of the transonic nozzle is 

A{x) = 1.398 + 0.347 tanh(0.8x — 4), x E (0, 10) (32) 

The inlet Mach number was set to M = 1.5, and the shock through the over- 
expanded nozzle occurred at x = 5.0. The resulting density profile from each 
ESWENO scheme is compared to the exact solution in Figure (5). The solution 
is non-oscillatory, and spreads the shock over a few cells. The results from the dif- 
ferent boundary treatments do not differ appreciably, which can be attributed to 
the zero gradient at the boundary. 


17 


1 


0.9 - 


0.8 


— Velocity 

— Pressure 
Density 

— Mach Number 



Figure 4. Solution for the quasi-ID subsonic nozzle problem is shown. 


Remark. The transonic one-dinrensional nozzle is a problem where boundary 
effects do not significantly influence the solution, but the problem demonstrates the 
stability of the ESWENO boundary closure for a steady state problem. The behavior 
of the solution to this problem in the current work exhibits the same characteristics 
that were observed in the periodic case [5] . Similarly, Sod’s shock tube problem and 
the sine-shock wave interaction problems exhibited the same characteristics. The 
results for those problems are not repeated here. 


5.3.3 Two Dimensional Inviscid Vortex 


An inviscid two-dimensional vortex convection problem was used to test the accu- 
racy of the finite-domain ESWENO framework in multiple dimensions. The exact 
solution to the convected vortex is 


f(x, y,t ) = 1 - ((x - x 0 - Uoct ) 2 + (y- ?/o) 2 ) , 

( o 7 — 1 \ 1 

,y,t) = ( 1 - e -^ 2 ~exp (f(x,y,t)) J , p = 

i .x jr y-yo ( f(x,y,t) 

x-xo-U^t ( f(x,y,t) 
v{x,y,t) = e exp' 


oi 

7 


(33) 


Uoo = MooCoo, x G (0, 10), y G (- 5,5), (x 0 ,y 0 ) = (5,0), t > 0, 

where e = 1.0, = 0.5, and 7 = 1.4 were used. Density contours from the initial 

condition and the approximate solution at t = 10 for the ESWEN 03 _ 4_3 operator 
with N = 200 x 200 grid cells are shown in Figure (6). 
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Figure 5. Density is plotted for Sod’s shock tube problem at t = 0.2 using N = 100 
cells. 


Both boundary closures were tested for the two-dimensional vortex and com- 
pared to corresponding linear operators. The accuracy is affected not only by the 
boundary closure of the finite difference scheme, but also by the treatment of the 
boundary conditions. In this case, the SAT penalty method is used to specify 
the exact boundary data in a well posed manner. This is meant to limit reflec- 
tions from the boundaries to be influenced only by the boundary closure and the 
SAT penalty term. A convergence study would not yield meaningful results using 
more common boundary condition treatments, such as Navier-Stokes Characteristic 
Boundary Conditions. 

The L 2 and convergence rates for the ESWENOg_ 4_3 operator are tabulated 
in Tables (6) and (7), respectively. The convergence rates for the ESWENC>2_4_2 
operator are given in Tables (8) and (9). The results show that the linear operators 
recover the design order accuracy. The ESWENO operators yield better than design 
order convergence rates for lower resolutions, but achieve design order convergence 
as resolution increases. The ESWENO operators exhibit larger errors for lower 
resolution than the corresponding linear operator. However, as the critical regions 
are better resolved, the linear operators and ESWENO operators yield similar L 0 0 
errors and nearly identical L 2 errors. This is the expected behavior, since the non- 
linear weights in the ESWENO operator approach the target weights as resolution 
is improved, recovering the linear operator. 


5.4 Supersonic Shear Layers 

A more practical test of ESWENO is the two dimensional non-reacting and react- 
ing supersonic hydrogen- air shear layers. These cases test the behavior of the new 
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(a) t = 0 


(b) t = 10 


Figure 6 . Density contours are shown for a two dimensional isentropic vortex solved 
using ESWEN 03 _ 4_3 


Table 6 . L 2 error and convergence rates for the ESWENC> 3 _ 4_3 and linear block 
norm centered difference operators for the two dimensional inviscid vortex. 



Linear Block Norm 

ESWENO 3 . 4.3 

Number of Cells 

L2 Error 

L2 Rate 

L2 Error 

L2 Rate 

50 x 50 

2.49E-05 

- 

5.32E-05 

- 

100 x 100 

1.64E-06 

3.92 

2.25E-06 

4.57 

200 x 200 

1.04E-07 

3.98 

1.08E-07 

4.37 

400 x 400 

6.57E-09 

3.99 

6.62E-09 

4.04 


scheme where shocks, high gradients, and physical instabilities are present and in- 
teract. A successful scheme must provide enough dissipation to damp oscillations 
near shocks and unresolved species and temperature gradients, but not be so dissi- 
pative that it damps the physical Kelvin-Helmholtz instability, which mixes the fuel 
and oxidizer streams. In addition to shocks and high gradients, the cases presented 
here include a discontinuous inlet condition and cannot be solved using standard 
centered finite differences without the addition of filtering and limiters. 

The hydrogen air mixing layer used for this study does not include any forcing 
of the inlet conditions. The top stream has molar species fractions Xh 2 = 0.5 and 
X 7 V 2 = 0.5, with a temperature of T = 300 IF. The bottom stream is composed of 
pure air, with inlet temperature of T = 300 K for non-reacting flow and T = 1500 K 
for reacting flow. The Mach number of the top stream is M = 2.0, while the bottom 
stream is M = 1.2. The full two dimensional compressible Navier-Stokes equations 
were solved. ESWENO was used to calculate convective gradients, while a 4th 
order, wide stencil SBP centered difference was used to calculate gradients of the 
viscous terms. The species were assumed to be thermally perfect ideal gases, and 
the open-source thermodynamics library Cantera [ 22 ] was used to calculate ther- 
modynamic and transport properties. A hydrogen air chemical kinetic mechanism 
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Table 7. L 0 0 error and convergence rates for the ESWENC> 3 _ 4_3 and linear block 
norm centered difference operators for the two dimensional inviscid vortex. 



Linear Block Norm 

ESWEN0 3 _4_ 3 

Number of Cells 

Loo Error 

Loo Rate 

Loo Error 

Loo Rate 

50 x 50 

2.46E-04 

- 

4.91E-04 

- 

100 x 100 

1.68E-05 

3.87 

3.31E-05 

3.89 

200 x 200 

1.07E-06 

3.97 

1.53E-06 

4.44 

400 x 400 

6.71E-08 

4.00 

9.04E-08 

4.08 


Table 8. L 2 error and convergence rates for the ESWENC>2_4_2 and linear diagonal 
norm centered difference operators for the two dimensional inviscid vortex. 



Linear Diagonal Norm 

ESWEN0 2 _4_2 

Number of Cells 

L2 Error 

L2 Rate 

L2 Error 

L2 Rate 

50 x 50 

3.49E-05 

- 

6.48E-05 

- 

100 x 100 

3.05E-06 

3.52 

4.43E-06 

3.87 

200 x 200 

3.40E-07 

3.17 

3.64E-07 

3.60 

400 x 400 

4.13E-08 

3.04 

4.16E-08 

3.13 


with 13 chemical species and 25 reactions was used. 

The density contours for the supersonic shear layers are shown in Figure (7). In 
both cases, the Kelvin-Helmholtz instability is clearly predicted, and the solution 
is stable despite the presence of shocks and unresolved species and temperature 
gradients. This suggests that the new ESWENO scheme exhibits the robustness 
required for chemically reacting calculations without adding excessive numerical 
dissipation. 

6 Conclusions 

A general strategy is presented in Refs. [5] and [6], for constructing Energy Stable 
Weighted Essentially Aon- Oscillatory (ESWENO) finite difference schemes on periodic 
domains. ESWENO schemes up to eighth-order are developed and proved to be sta- 
ble in the energy norm for both continuous and discontinuous solutions of systems 
of linear hyperbolic equations. Herein, boundary closures are developed for the 
fourth-order ESWENO scheme that maintain wherever possible the WENO stencil 
biasing properties, while satisfying the summation- by-parts (SBP) operator conven- 
tion, thereby ensuring stability in an L 2 norm. A novel set of non-uniform flux 
interpolation points is necessary near the boundaries to simultaneously achieve 1) 
accuracy, 2) the SBP convention, and 3) WENO stencil biasing mechanics. The 
novelty lies in the recognition that the discrete set of flux points 5cm nrust be consis- 
tent with the stability norm V used in the SBP formulation, and in fact are derived 
from the norm itself. Using the new flux-points fy, second-order, and third-order 
boundary closures are developed that achieve stability in a diagonal- and block-norm 
sense, respectively. 
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Table 9. Loo error and convergence rates for the ESWENC>2_4_2 and linear diagonal 
norm centered difference operators for the two dimensional inviscid vortex. 



Linear Diagonal Norm 

ESWEN0 2 _4_2 

Number of Cells 

Loo Error 

Lqo Rate 

Loo Error 

Loo Rate 

50 x 50 

3.17E-04 

- 

5.48E-04 

- 

100 x 100 

2.42E-05 

3.71 

4.46E-05 

3.62 

200 x 200 

2.28E-06 

3.41 

3.11E-06 

3.84 

400 x 400 

2.42E-07 

3.23 

2.74E-07 

3.51 


Extensive numerical validation is presented to assess the efficacy of the new 
boundary closures; the test problems included 1) unsteady 1-D linear wave equa- 
tion (hyperbolic), 2) steady quasi one-dinrensional nozzle flow solving the nonlin- 
ear Euler equations, 3) unsteady propagation of a 2-D Euler vortex, and 4) un- 
steady convection-diffusion-reaction of a supersonic hydrogen-air mixing layer. The 
global accuracy of the second-order diagonal-norm operators and third-order clo- 
sures block-norm operators, is shown to be “design-order”; i.e. three and four, 
respectively, for all test cases. Furthermore, it is shown that the accuracy of the 
new ESWENO operators is very close to that of the linear target operator. Finally, 
an eigenvalue solver is used to show that the ESWENO discrete operator is stable 
for smooth and discontinuous initial data. 

Four appendices are included to provide detailed instructions on the implemen- 
tation of the new ESWENO scheme. The appendices include sections describing the 

1) smoothness indicators used by ESWENO, including near- wall stencil definitions, 

2) a numerical recipe section that includes pseudo-code for the implementation of 
ESWENO on nonlinear hyperbolic wave equations, (i.e. Euler), 3) the diagonal- 
norm interpolation and weights formulae, and the dissipation coefficients needed to 
stabilize the baseline WENO approach, and 4) the block-norm interpolation and 
weights formulae, and the dissipation coefficients needed to stabilize the baseline 
WENO approach. 
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Density: 0.61 0.7 0.79 0.88 0.97 1.06 1.15 



(a) Non-Reacting Flow 



Density: 0.16 0.22 0.28 0.34 0.4 0.46 0.52 0.58 0.64 



(b) Reacting Flow 


Figure 7 . Density contours for a non-reacting and reacting supersonic shear layer 
utilizing fourth order ESWENO. 


Appendix A 

Smoothness Indicators 

~( r ) 

The smoothness indicators, , 2 , are calculated at all flux points using the 
expressions 

^i+1/2 = ( fj ~ fj- 1) Pf+1/2 = (/7+1 ~ fj ) ^j+1/2 = (fj + 2 ~ fj+ 1) 

1/2 = (/,-! - /,- 2) 2 ^ 1/2 = (/;+ 3 - / J+ 2) 2 

(Note that the indicators and j 3 RR are only used near the boundaries, and then 
only when a diagonal matrix norm is being used to prove stability.) 

To guarantee that the downwind stencil weight does not exceed that of the 
central or upwind weights, the downwind smoothness indicator is modified using 
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the expression 


/ \ 1 / fc 

^j+1/2 = ^3 [^j+1/2] j > ( A2 ) 

where k is an even integer. 

An additional stencil biasing parameter, fj+1/2, is needed for the ESWENO 
scheme. Here, fj+1/2 is a quadratic function of the highest undivided difference 
available on the stencil for fj+1/2 ■ From flux points Xf 2 to Xf N2 [see Fig. ( 1 )], 

fj+ 1/2 = (~fj - 1 + 3 fj — Sfj+i + fj+2 ) 2 (A 3 ) 

At the point Xf t , fj+1/2 is biased toward the interior of the domain, 

fj+ 1/2 = (~fj + 3 /j+i — ^ fj+2 + fj+3 ) 2 j (A 4 ) 

and at Xf N _ 1 , fj+1/2 is biased in a mirrored fashion. 

The parameter e is a function of number of points in the discretization. 

e = max(\\fo\\,\\fo\\) x ^ Xd (A£) 3 , (A 5 ) 

where ||/o|| and || IK represent a norm of the flux and the gradient of the flux from 
the initial condition, where points near a discontinuity are excluded. A1 


1 For applications where the solution changes drastically during simulation, e can be rescaled as 

necessary. 
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Appendix B 


Recipe 


The recipe for calculating the gradient using the fourth-order ESWENO schemes 
is summarized below. 

1. Construct the convective flux f at solution points (xjy) according to the gov- 
erning differential equation(s). 

2. For systems of equations (e.g. Euler), construct the flux Jacobian matrices: 
A = gjj at the flux points ( xm ) using Roe average variables formed from 
nearest neighbor solution point data. Form the eigenvector decomposition 

A = SAS- 1 . 

3. At the flux point x-,i, identify all solution points x^r that contribute to the 

2 

r r \ 

interpolants , ?' = L, C , R. Use the eigenvector matrix S . j to rotate the 

•7 '3 

solution and flux at these points, into characteristic form C = S~ l U ; F c (C ) = 

s^Fiu). 

4. Form the Lax-Friedrichs characteristic fluxes fA = | (f c ± A maa; c). 

5. Perform interpolations on each candidate stencil, = mT\ f c . 

6. Calculate the stencil biasing parameters: 

(a) Calculate /3 r for each flux point according to Eq. (Al), except at the 
endpoints. 

(b) Calculate fj+ 1/2 according to Eq. (A3) and Eq. (A4). 

7. Calculate and normalize the weights using the stencil biasing parameters and 
the target weights in Eq. (14). 

8. Calculate A* from the weights using Eqs. (C3), (C4), and (C5) for the diagonal- 
norm scheme or Eqs. (D5) and (D6) for the block-norm. 

9. Modify the diagonal of Aj to be smoothly positive according to Eq. (20). 

10. Calculate the WENO flux from the weights and candidate interpolations using 
Eq. (13). 

11. Calculate the Energy Stable flux, if;, from Eq. (23). 

12. Reconstruct the fluxes in characteristic space, fj+ 1/2 = fj+ 1/2 + fj+ 1/2 anc ^ 
V^j+ 1/2 = V^j+ 1/2 + V’j+ 1 / 2 i an d transform back to physical space. 

13. Calculate the gradient using the inverse of the P-norm, as in Eq. (24). 
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Remark. It is important to note that all equations presented are for the forward 
propagating waves, The equations for interpolation of backward propagating 

waves (,/j+i / 2 ) are found simply by flipping r = (LL, L, C, R, RR ) about C. For the 
A i terms, after flipping r about (7, the diagonal matrix is multiplied by —1. Eq. 
(23) is also multiplied by —1 to calculate 
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Appendix C 


Fourth-order Diagonal-Norm Operator: D 2-4-2 


C.l Differentiation Matrix 


The target SBP differentiation operator T>2-4-2 


D = - — 


1 

Ax 


( 

24 

I, 7 

59 

34 

4 

47 

3 

34 

0 

0 

0 

0 

0 

0 

0 
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2 

0 

1 

2 
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0 

0 

0 

0 

0 

0 

0 


4 

43 

59 

86 

0 

59 

86 

4 

43 

0 

0 

0 

0 

0 

0 


3 

98 

0 

59 

98 

0 

32 

49 

4 

49 

0 

0 

0 

0 

0 


0 

0 

1 

12 

2 

3 

0 

2 

3 

1 

12 

0 

0 

0 

0 


0 

0 

0 
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0 
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0 

1 

12 

2 

3 
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12 
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0 

4 

49 

32 

49 

0 

59 

98 

0 

3 

98 


0 

0 

0 

0 

0 

0 

4 

43 

59 

86 

0 

59 

86 

4 

l 43 


0 

0 

0 

0 

0 

0 

0 

0 

1 

2 

0 

1 

2 

V 

0 

0 

0 

0 

0 

0 

0 

3 

34 

4 

17 

59 

34 

24 

17 


where the diagonal mass matrix V is 


P = AxDiag (I | I | 1 


i 49 43 59 17 \ 

1 48 48 48 48 ) 


The skew-symmetric matrix “Q” follows immediately from P D. 
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C.2 Interpolation Operators 


The WEN0/ESWEN02_4_2 requires the inclusion of two extra stencils near the 
boundary. The interpolation coefficients are presented below. 



C.3 Target Weights 

The target weights for the five interpolation stencils are: 


d-242 



0 

1 

0 

0 

11 

24 

1013 

31 

51 

48^8 

& 

70 

408 

? 

32^66 

6 

5p 

f 

6 

3 

1 

2 

1 

$ 

4§8 

3§7 

f 

5 5 7 1 5 

32(j>6 

10?3 

52 

56 

0 

4898 

31 

0 

1 

0 



(C2) 
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C.4 Energy Stable Terms 


The nonzero terms of the diagonal matrices, A,;, for the ESWENC>2_4_2 are 


/ 


0 

± (l4 - 28 w { 2 L) - 71 wf L) - 31*3 - 3 lw[ RR) ^j 
± ^—8 + 28 w ( 2 L) - 23u4 L) + 48w ( 3 LL) + 31*3 - 20*3 - 48 w[ RR) ^ 
^{2 + 23 - 24w{ L) + 23 wf L) + 20iv (R) - 2bw {R) + 79 w[ RR) ^j 
± (24w{ L) - 24 wf ] + 25w { 3 R) - 24 w[ R) ^j 


Diag(Ai) = 


1 (..& _ W ( L ) 

i w i + 1 


W 


(H) „.,(«) 


i— 1 


— it;; 


96 ( 24u, AT-4 25 3-3 + 24u, AT-5 24u, AT-4) 

gg (-2 + 25u>3 - 79*43 - 20*43 + 24 *43 " 23*43 " 23*43) 
<k ( 8 - 31*43 + 48*43 + 20*43 + 23*43 - 48*43 - 28*43) 
M (" 14 + 31*43 + 3l4S + 71 w {RR l + 28 *43) 

V ° 

(C3) 


( 0 \ 

g ' 6 (28*4 L) + 142*4 LL) - 31*4 R) - 110*4 RR) ) 

^ (3*3 + 94u 3 - 20u;^ ) - 158^3) 
m ( 24w 4 L) - 25*3) 


Diag (A 2 ) = 


1 

4 



— w 


(R) 
1- 1 


(C4) 


— f 25**/^ — 24 ^ 

96 \^ OW N-3 ^ W N- AJ 

<k (158^3 + 20*43 - 23*43 - 94**3) 

<k (3l3i + 1 10*43 - 142 **3 " 28*43) 

V 0 / 


29 



/ 


96 



0 

0 

(LL) 


+ 79 w[ RR) ^j 


Diag (A 3 ) = 


0 


96 


V 


(— 79 ^ 71 ^) 

0 

0 


The matrix Diag (Aq) is identically zero. 
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Appendix D 


Fourth-order Block-Norm Operator: F ) 3_4_3 


D.l Differentiation Matrix 


The target SBP differentiation operator satisfies the expression T ) 3 _ 4_3 = 

^3-4-3^3-4-3 [ see ec L ua ti° n (10)]- The inverse P-norrn for the ESWENC> 3 _ 4_3 scheme 
is given by 


P. 


-1 

3-4-3 


1 

Ax 


/ Po 1 

0 

0 \ 


0 

I 

0 

p~ 1 — 
> M) ~ 

V 0 

0 

(p 0 -r T ) 
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m 4 85 4 3 
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1 
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where PT denotes the per-symmetric transpose. 


V 

(Dl) 

The Q 3 _ 4_3 matrix is given by 
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(D2) 
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D.2 Interpolation Operators 


The interpolation coefficients for the ESWENC>3_4_3 are 
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(D 3 ) 


D.3 Target Weights 


The target weights for each stencil in ESWEN03_4_3 are given below for each flux 
point. 
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D.4 Energy Stable Terms 

The nonzero terms of the diagonal matrices, A,;, for the ESWENC>3_4_3 are 


Diag(Ai) = 


2gg (58 - 100'«4 L) - 101m} R) ) 

2§g (—56 + lOOm^ 4 — Glw^ + lOlm^ — 44^4^^ 
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